/* 
 * File:   main.cpp
 * Author: jcfernandez
 *
 * Created on 21 de noviembre de 2011, 16:25
 */

//******************************************************************
//  MDDipFlow.cpp
//  Programa principal de simulacion de particulas con momento
//	magnetico constante bajo un campo magnetico externo
//  El sistema se encuentra confinado entre paredes situadas en el eje z
//  en posiciones +Lz/2,-Lz/2
//
//	Grupo de Fisica de Fluidos y Biocoloides - Univers. de Granada
//  �2011, Dr. Juan Carlos Fernandez Toledano
//******************************************************************
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>
#include<iomanip>
#include <sstream>

//#define CALC_STRESS
#define POLYDISPERSE
#define CONFINED
//#define CORRELADORES
#define PRINTXYZ
//#define FLUJOS

using namespace std;

#define kB 1.0//8.314e-03
#define pi 3.1415926535898
#define avogadro 602.205

template<typename T>inline T SQR(T x) {
    return (x * x);
}

template<typename T>inline T POW3(T x) {
    return (x * x * x);
}

template<typename T>inline T VLENGTH3(const T *x) {
    return (sqrt(SQR(x[0]) + SQR(x[1]) + SQR(x[2])));
}

template<typename T>inline int anint(T x) {
    return (int(floor(x + 0.5)));
}

#include "utils.h"
#include "particulas.h"
#include "caja.h"
#include "correlador.h"
#include "interaccion.h"
#include "verletClass.h"
#include "verletBrownian.h"
#include "espreso.h"
#include "print.h"

int main(int argc, char** argv) {

    ///Tiempo de la simulacion
    double t;
    ///Parametro gamma de la resolucion de las ecuaciones de movimiento
    double gamma = 10.0;
    ///Empleo del termostato
    bool termo = false;
    ///Caja de la simulacion
    caja *box;
    ///Correlador
    correlador *correlator;
    ///Imprimimos resumen simulacion
    printTilte();
    ///Cargamos el archivo de configuracion
    parametros param;
    param.importConfig();
    param.printConfig();
    ///Nombre del archivo de tensor de esfuerzo
    char name[50];
    ///Cargamos el nombre del fichero del tensor de esfuerzos
    strcpy(name, param.stressFile);
    ///Lo borramos si existe
    ofstream fres(name, ios::out);
    fres.close();
    ///Semilla aleatoria
    //srand (time(NULL));



    ///Creamos la caja de simulacion
    box = new caja(param);
    ///Creamos el archivo binario de almacenamiento
    ofstream ffout("posvel_bin.dat", ios::out | ios::binary);
    ffout.write((char *) (&box->npart), sizeof (int));
    ffout.write((char *) (&box->lbox[0]), sizeof (double));
    ffout.write((char *) (&box->lbox[1]), sizeof (double));
    ffout.write((char *) (&box->lbox[2]), sizeof (double));
    ffout.close();

#ifdef CORRELADORES
    ///Creamos los correladores
    correlator = new correlador(box, param);
    ///Inicializamos el fichero de salida de los correladores
    correlator->initFileCorr();
#endif

    //verlet *vrl= new verlet(box,Dt,termo);
    //verlet *vrl= new verlet(box, param.Dt, gamma, termo);
    double temp = 1.0;
    //verlet *vrl = new verlet(box, param.Dt, temp, 0.05);
    //verletBrownian *vrl = new verletBrownian(box, param.Dt, temp, 0.005);
    verletEspresso *vrl = new verletEspresso(box, 10.0, 2.0, param.Dt);

#ifdef CONFINED
    ///Calculamos la interaccion del sistema con confinamiento
    interaccionConfined(box);
#endif

#ifndef CONFINED
    ///Calculamos la interaccion del sistema sin confinamiento
    interaccionSimple(box);
#endif

    ///Creamos el archivo de energia
    ofstream fout("energia.dat", ios::out);
    fout.precision(5);
    fout.is_open();


    ///Comienza el bucle principal
    t = 0.0;
    for (int i = 0; i < 100000000; i++) {

        ///Actualizamos el paso de tiempo
        t += param.Dt;
        box->nacc = i;
        //vrl->PaulYoonEngine();
        //vrl->verletLangevin();
        vrl->integrate_vv();

#ifdef CORRELADORES
        correlator->calcCorrelador(t);
#endif

        if (i % 1000 == 0) {

            //time_t tAct = time(NULL);
            printf("%10.5f  %8.2f  %8.2f  %8.2f  %8.2f\n", t, box->epoten, box->ekin, box->epoten + box->ekin, box->temp);
            fout << t << "  " << box->epoten << "  " << box->ekin << "  " << box->epoten + box->ekin << "  " << box->temp << endl;
            //printConfigBackup(box,param);
            box->printConfigBackup();

#ifdef PRINTXYZ
            box->printXYZfile(param.confFile);
            //printXYZfile(box);
#endif

#ifdef CALC_STRESS
            box->printShearStressTensor();
            //printShearStressTensor(box, param);
#endif
        }
    }
    fout.close();
    delete []box;
    delete []correlator;


    cout << "sacabo" << endl;
    return 0;
}
